# -*-coding:utf-8 -*-
# @Time: 2023/3/27 16:03
# @Author: cuishuohao
# @File: demo
# @Software: PyCharm

# 1.构建多分类数据集
# (1) 首先生成4个分类的数据
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(0)
Num = 100

# 第1个分类的数据
x_1 = np.random.normal(-3, 1, size=Num)
x_2 = np.random.normal(-3, 1, size=Num)
y_0 = np.zeros(Num)  # 这里的y指的是该样本属于正例的概率，全为0
c_0 = np.array([x_1, x_2, y_0])  # 0分类

# 第2个分类的数据
x_1 = np.random.normal(3, 1, size=Num)
x_2 = np.random.normal(-3, 1, size=Num)
y_1 = np.ones(Num)
c_1 = np.array([x_1, x_2, y_1])

# 第3个分类的数据
x_1 = np.random.normal(-3, 1, size=Num)
x_2 = np.random.normal(3, 1, size=Num)
y_2 = np.ones(Num) * 2
c_2 = np.array([x_1, x_2, y_2])

# 第4个分类的数据
x_1 = np.random.normal(3, 1, size=Num)
x_2 = np.random.normal(3, 1, size=Num)
y_3 = np.ones(Num) * 3
c_3 = np.array([x_1, x_2, y_3])

# 调整顺序
c_0 = c_0.T
c_1 = c_1.T
c_2 = c_2.T
c_3 = c_3.T

# 图形显示
plt.scatter(c_0[:, 0], c_0[:, 1], c='b', marker='x')
plt.scatter(c_1[:, 0], c_1[:, 1], c='r', marker='+')
plt.scatter(c_2[:, 0], c_2[:, 1], c='y', marker='v')
plt.scatter(c_3[:, 0], c_3[:, 1], c='g', marker='s')
plt.xlabel("x_1")
plt.ylabel("x_2")
plt.show()

# (2) 构建训练集和测试集
All_data = np.concatenate((c_0, c_1, c_2, c_3))
np.random.shuffle(All_data)  # 打乱数据
# 分别构建训练集和测试集
train_data_X = All_data[:300, :2]
train_data_y = All_data[:300, -1].reshape(-1, 1)
test_data_X = All_data[300:, :2]
test_data_y = All_data[300:, -1].reshape(-1, 1)

# 2.建立Softmax回归模型
# (1)由于有4个分类，因此需要初始化4个模型的参数
# 初始化W和b，为了方便展示，注意这里没有将参数初始化为0
# 4行2列的随机数
w = np.random.rand(4, 2)  # 4个线性函数，特征x1和x2各需要一个参数w
# 1行4列的随机数
bias = np.random.rand(1, 4)  # 4个线性函数各自需要一个 bias

# 展示初始的分类器结果:
plt.scatter(c_0[:, 0], c_0[:, 1], c='b', marker='x')
plt.scatter(c_1[:, 0], c_1[:, 1], c='r', marker='+')
plt.scatter(c_2[:, 0], c_2[:, 1], c='y', marker='v')
plt.scatter(c_3[:, 0], c_3[:, 1], c='g', marker='s')
plt.xlabel("x_1")
plt.ylabel("x_2")

# 含左不含右，步长为1的 ndarray
x = np.arange(-6, 6)
y1 = -(w[0, 0] * x + bias[0, 0]) / w[0, 1]
plt.plot(x, y1, "b")
y2 = -(w[1, 0] * x + bias[0, 1]) / w[1, 1]
plt.plot(x, y2, 'r')
y3 = -(w[2, 0] * x + bias[0, 2]) / w[2, 1]
plt.plot(x, y3, 'y')
y4 = -(w[3, 0] * x + bias[0, 3]) / w[3, 1]
plt.plot(x, y4, 'g')
plt.show()


# (2)定义相关函数
# 定义softmax函数，如果输入是向量则有问题，应该如何改进 softmax函数?
def softmax(z):
    sum_exp = np.sum(np.exp(z), axis=1, keepdims=True)
    return np.exp(z) / sum_exp


# 对train_data_y 进行one-hot编码
def one_hot(temp):
    one_hot = np.zeros((len(temp), len(np.unique(temp))))
    one_hot[np.arange(len(temp)), temp.astype(np.int16).T] = 1  # 注意这里的y需转换为整数，否则会报错
    return one_hot


# 计算线性模型结果
def compute_y_hat(W, X, b):
    y_hat = np.dot(X, W.T) + b
    return y_hat


# 根据公式3.44定义损失的数
def cross_entropy(y, y_hat):
    return -(1 / len(y)) * np.sum(y * np.log(y_hat))


# (3)使用训练集训练模型
loss_list = []
epochs = 1000
lr = 0.01
# 训练模型
for i in range(epochs):
    y = one_hot(train_data_y)
    y_hat = softmax(compute_y_hat(w, train_data_X, bias))
    loss = cross_entropy(y, y_hat)
    loss_list.append(loss)
    # 计算梯度
    # grad_w: 2行4列
    grad_w = -(1 / len(y)) * np.dot(train_data_X.T, (y - y_hat))
    # grad_b: 4个元素
    grad_b = -(1 / len(y)) * np.sum((y - y_hat), axis=0)
    # 更新w和b
    # w: 4行2列
    w = w - lr * grad_w.T
    # bias: 1行4列
    bias = bias - lr * grad_b
    if i % 100 == 0:
        print("i:", i, "  loss:", loss)
    if loss < 0.0001:
        break

# 查看 Loss值的变化情况:
plt.plot(loss_list)
plt.show()


# (4)查看模型在测试集上的表现
def predict(x):
    y_hat = softmax(compute_y_hat(w, x, bias))
    return np.argmax(y_hat, axis=1)[:, np.newaxis]


# 输出结果是 1.0 说明该模型可以正确分类所有的测试集数据。
print(np.sum(predict(test_data_X) == test_data_y) / len(test_data_y))
